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Abstract: We examine both analytically and numerically the validity of factorization for 
the double dipole scattering amplitude T^'^^ which appears on the right hand side of the 
BK-JIMWLK equation. We demonstrate that, if one uses a dilute object (e.g., a proton in 
DIS) as the initial condition, the correlation in the transverse plane induced by the leading 
order BFKL evolution is generally strong, resulting in a violation of the mean field approx- 
imation T*^^) ~ TT even at zero impact parameter by a factor ranging from 1.5 to 0(10) 
depending on the relative size of the scatterers and rapidity. This suggests that, within 
the experimentally accessible energy interval, the transverse correlation can significantly 
aff'ect the nonlinear evolution of the dipole scattering amplitude. It also suggests that the 
nonlinear eff'ects may set in earlier, already in the weak scattering regime. In the case of 
the simulation with a running coupling, the violation of factorization is somewhat milder, 
but still noticeable. 
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1. Introduction 

High energy scattering near the unitarity limit is a delicate problem which deserves intense 
theoretical efforts in view of its phenomenological importance at hadron colliders. There is 
a clear goal of including nonlinear, saturation effects due to the high density of gluons into 
the energy evolution of scattering amplitudes, but a precise determination of when and 
how these effects should be treated is subject to various uncertainties depending on the 
process of interest. The problem appears to somewhat simplify if one considers scattering 
of a small object (e.g., a photon at high virtuality in DIS) off a very heavy nucleus where 
saturation is important already at relatively low energy. For such a process the Balitsky- 
Kovchegov (BK) equation [1,2] is the most commonly studied equation which provides a 
concrete scenario for an approach towards unitarity, 

ayTy(x,y) = ^ y d2zX(x,y,z)|-ry(x,y)+ry(x,z) + ry(z,y)-ry(x,z)ry(z,y)|, 

(x-z)^(.-,)^ ' «^ = ^- (1-1) 

Here Tyix^y) is the forward amplitude of a dipole of size |x — y| at rapidity Y . The first 
three terms on the right hand side contain the BFKL physics [3, 4] while the last term 
~ TT ensures that the amplitude saturates the black disc limit T — > 1 which is a fixed 
point of the equation. Being a closed equation, ( |1 . 1[) is amenable to both analytical and 
numerical approaches, and the properties of the solution as well as their phenomenological 
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consequences have been discussed extensively over the past several years (see, reviews [5,6] 
and references therein). 

However, it is not often emphasized that the BK equation is a mean field approximation 
to a more general equation, namely, the B-JIMWLK equation [1,7-10] 



dvTYix.y) = — J d zM{x,y,z)^-TY{x,y) +Ty{x,z) +TY{z,y) - {Ty {x , z)Ty {z , y)) 

(1.2) 

nor is the validity of this approximation fully appreciated. Here the brackets (...) denote 
averaging over the target configurations. The difference between these two equations is 
usually considered to be minor: Although the former obviously discards any kind of ex- 
isting correlations in the target wavefunction, this would be justified for a large nucleus 
at low rapidity (see, however, [11]). The subsequent quantum evolution then generates 
correlations which vanish in the large Nc limit, 

{TT)^{T){T) + o(^-^^ . (1.3) 

Indeed, the only existing numerical simulation of the B-JIMWLK equation [12] starting 
from uncorrelated initial conditions shows little difference from the corresponding solution 
to the BK equation. 

The purpose of this work is to demonstrate that the factorization ( |1.3| ) is violated 
when one considers a dilute target consisting of a few partons (e.g., a proton) instead of a 
heavy nucleus as the initial condition. Of course, there is a priori no reason to expect that 
factorization should work in this case, but there has not been any quantitative study of the 
degree of its violation either. For a dilute target, a significant part of the rapidity evolution 
in realistic experiments is in the linear BFKL regime where the amplitude is rapidly growing 
but still much less than unity, whereas saturation is considered to be relevant only in the late 
stages of the evolution.^ The fluctuations and correlations developed in the linear regime 
are so strong that the initial condition that should be used for the nonlinear evolution 
equations is a highly nontrivial system of gluons for which the difference between ( |1.1| ) 
and ( |1.2| ) may turn out to be crucial, especially for phenomenology. Specifically, in the 
framework of the QCD dipole model ref. [13] found a power-law correlation in the double 
scattering amplitude^ 

{T{x,z)T{w,y)) ^j^^^, (1.4) 

under the condition that the distance between the two dipoles are much larger than their 
sizes, |z — wl » \x — z\, \w — y\. (7 is a positive, calculable number related to the 
anomalous dimension.) In the exemplary cases studied in [13], this power-law always leads 



^However, we have found some evidence that nonlinear effects might set in earlier due to the correlation. 
See the discussion in section 



3.2 



"See also [14], though there seem to be disagreements in the results. 
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to a parametrically large ratio 

{T(x,z)){T(w,v)) ^ ' 

Due to a technical reason, in [13] it was not possible to take the interesting limit w ^ z 
to evaluate R for the 'BK configuration', although it was tantalizing to conclude from ( |1.^ ) 
that that the correlation would become even larger in this case. Here we circumvent this 
difficulty and present an analytical insight into the behavior of as a function of the initial 
dipole sizes. 

However, analytical calculations are often quite difficult, and one can usually only deal 
with special configurations which are set by hand. Besides, for our purpose it is important 
to know the actual numerical value of T and (T^) to make sure that one evaluates R in 
a regime where the nonlinear corrections just start to be important. We will therefore 
also perform a Monte Carlo (MC) simulation of the QCD dipole model [15] which contains 
the exact leading order BFKL dynamics. In this framework one generates a cascade of 
dipoles keeping track of their sizes and positions in the transverse plane. Calculations of 
{T^) for any k, hence R, are completely straightforward for arbitrary configurations. We 
then compare the numerical results with analytic expectations and find that they agree 
satisfactorily. For zero impact parameter we find that R is much larger than 1 when the 
ratio of the projectile and target sizes is either small or large. The minimum value for R is 
attained when the projectile and target are of similar size, and in this case the value of R 
is around 1.5. This suggests that, in the leading logarithmic approximation on which both 
the BK equation and the dipole model are based, the replacement (TT) (T)^ is not 
valid for a proton target especially for a small dipole projectile (or in the high-Q^ region 
of DIS), although it might be safe to do so for a nucleus target. In the former case one 
should rather use the B-JIMWLK equation with a strongly correlated initial condition, 
whose asymptotic solution can be different from that of the BK equation. 

The fact that one finds large correlations in the leading order evolution for a dilute 
system is consistent with the early studies on fiuctuations in [16, 17]. In [16] it was found 
that (r'^) ~ (^0^ (or rather (T^) ~ A:! • (A: + 3)!) at zero impact parameter. This implies 
that, for any m < k, 

{T'') ( kl Y (k^^ 



»1. (1.6) 

Note, however, that the definition of (T'^) in ( |1.6| ) is different from the one considered in 
this paper, namely, {T^) appearing in the Balitsky hierarchy whose first equation is (|1.2|). 
In (^), one evolves the target and the projectile up to some energy, and then calculate the 
sum of all events in which there are k simultaneous interactions. In our case we rather fix 
k given dipoles in the transverse plane, and then consider their scattering off some target. 
Only the latter contains information of the correlation resolved in the transverse plane. 

In [18-20] the dipole model has been modified and extended to include various non- 
leading effects as well as saturation and confinement effects during the evolution. Generally 
speaking, these effects tend to reduce the correlation. For example, {T'^) as defined in [16] 
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behaves as (for k between 5 and 9) {T^) / {T^-'^) ^ 1.2 ■ k once the nonleading effects are 
included [20]. This imphes 



(1.7) 



and thus the correlation is reduced with respect to (1.6). It should, however, be said that 
the fluctuations are still very important, and they have for example important consequences 
on the study of elastic and diffractive scattering in DIS and pp collisions [20]. In this paper 
we only show some of the preliminary numerical results with the running coupling effect to 
see if there is a similar suppression of the correlation, while a detailed study of the various 
additional effects is postponed to a future publication. 

The paper is organized as follows. In the next section we present analytical calcula- 
tions of the double dipole scattering amplitude and the ratio R for the BK configurations 
mentioned above. In section p.l| we outline our numerical approach to the calculation of 



the correlation. The results, including the running coupling case, are then presented in 



section 3.2 where we also make comparison with the analytical expectations. Finally, in 



section Q we summarize our results and raise some open questions. 

2. Analytical approach 
2.1 The dipole pair density 

In the dipole model [15], the degree of the two-body correlation in impact parameter space 
is encoded in the dipole pair density [21,22] whose integral representation reads (keeping 
only the zero conformal spin sector) [23, 24] 



d^X2d^X'id^Xi^ Ai7^1-7a/ \Tpl-lb( \ 
2 2 2 ^ ^{X2j,X3.y)E ^''{X2a,X4a)E 2^3/3 , 2:4/3) 



2'232'342'42 



(2.1) 



where xqi = xq — x\ denotes the coordinate of the parent dipole, and Xa^ax = Xa^ — Xa^ 
and Xfegb^ = Xfeg — x^^ are those of the child dipoles (see. Fig. ^. We shall use the letter x 
for both two-dimensional real vectors and their magnitude, x is the BFKL eigenvalue 

Xici) = a,(2V(l) - iPi-f) - ^(1 - 7)) , (2.2) 
with 7 being the anomalous dimension, and E is the eigenfunction of the SL(2,C) group 

E^{xo.y,xi^) = (^^y'' . (2.3) 
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Figure 1: A graphical representation of equation (2.1) 



The 7-integrals are along the imaginary axis. With the usual representation 7 = ^ + iu, 
it reads 



(2.4) 



In ref. [13], the multi-dimensional integral in ( |2.1| ) has been carried out in the limit 



^ab — ^b 



^ Xagai ; Xbgbi 



2.5) 



The result shows a power-law correlation between the two child dipoles. In the case of 
xoi ^ Xab, ref. [13] found 



^ab 



{ny 



(2.6) 



where n is the single dipole density, and ja and 7 are the saddle point values determined 
from certain conditions. The breakdown of factorization is carried over to that of the 
two-dipole scattering amplitude 



' \ 2{27a-7) 

T^'^\xaoai,Xbobi) ~ ( T{Xaoai)T{Xbobi) » T{Xaoai)T{xh^bi) , 



Xab J 



(2.7) 



as already noted in the introduction. (From now on we use the notation T^^^ in place 
of (T^).) On the other hand, the quantity of interest for us is the two dipole scattering 
amplitude for contiguous dipoles, namely. 



Xa\ — Xb-^ 



(2. 
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Although it is not legitimate to extrapolate the result (|2.q ) to the ^ 0, it does 

suggest that the correlations would become even larger for such 'BK configurations'. (The 
numerical evaluation of this case is presented in section ^^.) In this section we attempt 
at an analytical evaluation of n^^) for Xa^ = Xb-^ in certain limits and discuss the behavior 



of the ratio R defined in (1.5). The result will be confronted with numerical Monte Carlo 



simulations in the next section. 

2.2 Calculation of n^^^ for contiguous dipoles 

The last line of (^]^) is a known integral whose overall structure is fixed by conformal 
symmetry. After performing this integral, the last two lines of (^]^) become 



I = f{l,la,lb) / d Xad X/3 



d2^ ( \ \ (2 9) 

where the function /-the 'triple Pomeron vertex'- can be found in [25,26], and we have 
already set Xa^ = Xbi = Xc- 

To make progress we assume that 7^ = 76, which is a good approximation when the 
configuration of the two child dipoles is more or less symmetric. (The saddle points 7^ and 
7b depend only logarithmically on dipole sizes.) Then the x^ integral can be done [27] 

27 / ^ „ \ 2(1-7) 



+ (7 ^ 1 - 7)V (2.10) 



2(1-7) / ^^^■^ 



where F is the hypergeometric function. 



vr2-^-i4i±^^:Hz!^, (2.11) 



^2 — iy)T{l + iv 
and 

(2.12) 

is the anharmonic ratio of the four points (xq, xi, x^, x^j) {z is the complex coordinate 
representation of x), see fig. ||. The remaining integrals are difficult to perform in full 
generality. As in [13], we shall restrict ourselves to two limiting cases xqi — >• (small 
parents) and xqi — >■ 00 (large parents). In both limits, \p\ <^ 1, so we may approximate 
F(...,p) ~ 1. The two terms in ( |2.10| ) give equal contributions due to the symmetry 
7^1—7. Taking this into account, we can write 

/ = 2c,/(7,7a,7a) / ( ^^V'° ( Vf^V'" ( ^^Y" . (2.13) 

J Xli^ \Xa0aXci3J \Xb0i3XcaJ \XOaXipJ 
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The integrand is a product of anharmonic ratios weighted by the conformally invariant 
measure (Pxad'^xp/x^^, so it is invariant under conformal transformations of the external 
points. However, since there are five of them {xq, xi, Xa^, x^g and Xc), conformal symmetry 
is not strong enough to constrain the solution, and our assumption xqi — > oo or xqi — > 
will be crucial in the following. 



2.2.1 Large parents 

Suppose the parent dipole is large and the points XaQ,bo,c are all located near the center of 
the parent dipole as illustrated in fig. §(a). This may be regarded as a situation relevant to 
DIS on a hadron at high photon virtuality. Without loss of generality, we can set Xc = 0. 
The integrand vanishes as Xa,i3 — ^ oo very fast, so that a finite region of Xa,i3 near the origin 
is important. Therefore we may approximate 



xoi 4 



X0aXli3 xoi 
Under this assumption, ( ^.13 ) takes the form 



(2.14) 



\XmJ J xli^^ \XaoaXpJ \Xb^(}XaJ 

For simplicity, we assume that the two dipoles have the same size: Xa^c = Xb(^c = f- (The 
region Xa^c ^ Xb^c or Xq^c Xb^c gives a subleading contribution in the BFKL or the BK 
equation, see Section 2.3.) Writing z^^ — re* " and Zbg = re*^'' and rescahng x^,/? rxa^p 
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Figure 3: Graphical representation of the 'BK configurations', (a) for the large parent, small 
impact parameter case, and (6) for the small parent, large impact parameter case. 



we get 

= 2c^/(7,7„,7„)(— ) \{e), (2.16) 

where 6 = 9a — 9fj is the relative angle between the two child dipolcs. Wc have not been 
able to determine the function g(9) for ^ 7^ in a closed form (^(O) is a known integral 
in the conformal field theory literature [28-30]). But since this function has no singularity 
and depends only on the angle, it will not affect the evaluation of the saddle point below. 

Neglecting this angular dependence and other prefactors, we can estimate the two 
dipole scattering amplitude as 



^ {xqi, XaQCT XjjQc) ~ '^s-^a()C'^6()c'^'' (-^Ol j ^^aoC) ^^^yc) 

\xoiJ 



~ 1^ f djd-fad-fb ( ^ e^{'r)y+(x(ia)+x(it)Ky -y) . (2.17) 



'0 

After performing the y integral, we get the two contributions 

27 



/ d7d7„d7b — — , (2.18) 



and 

27 Jx{ia)+xht))y 



r / r \ e(xha}+xhb))y 

J \X0lJ Xila) + Xilb) - XiV 
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The saddle point for the 7a and 7f, integrals in ( 2.19| ) is simply the BFKL one 7a = 76 = 1/2, 
leading to 



where 7 solves 



r \ 



27 p2x(l/2)F 



XQiJ 2x(l/2)-x(7) \xoiJ 



27 



=2x(i/2)y 



X(7) = 2x 



7 0.82. 



For the contribution ( 2.1g| ) we can use the saddle point for the 7 integral, 



(2.20) 



(2.21) 



and the leading rapidity behavior of this contribution is then given by 



\xoiJ 



\ 27s 



(2.22) 



(2.23) 



As we discuss in section 2^, it holds that 2^(1/2) > xils), ^-e., 7s < 0.82 for all config- 
urations we are interested in. (In the limit Y — s- 00, 7^ 1/2.) The contribution which 
dominates is thus given by (|]20|), and we therefore have 



». \27 
!_ ] g2x(l/2)y 



xoij 



On the other hand, the single dipole scattering amplitude is given by 

27 



r(xoi,r)~a, f— ) e>^(^)^ 



\xoiJ 



where 7 is the solution to 



(2.24) 



(2.25) 



X'(7)l' = ln^. 



(2.26) 



Taking the ratio, we arrive at 



R=l!-L^ ^X01^, 2(27-7) 2(x(l/2)-x{7))> 
— — \2 



(T)2 



\ r J 



(2.27) 



Since 27 > 1 > 7, the first factor is larger than 1 and predicts that the correlation 
increases as the asymmetry becomes larger xqi » r. Since xil) > x(l/2), the second, 
exponential factor tends to decrease the correlation at high values of rapidity. Comparing 
this with ( |2.7| ), we infer that R monotonously increases and eventually saturates to the 
expression (|2.27| ) as Xab 0. 
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2.2.2 Small parents 

Another tractable example is the limit of a small parent dipole xqi — > 0. In this case 
we may approximate xip ~ XQj^, after which the point xi drops out from the integral. 
Rewriting 

r n J!/ . f XoiXaobo\'^'^ f(fxa(fxp( XaocXbgcXlf^ \ '^'^ ( XapXQaoXQbo^'^'^ 



\xQai^XQf)fj J J "^afl \ XaQaXcaXbopXclS I V ''^Oo'^O/J'^aotio 



(2.28) 



we see that, apart from the prefactor, the integrand is conformally invariant, so it can be 
written as 

/ = 2c,/(7, la, la) f h (^, ry*) , (2.29) 

where r] is an anharmonic ratio 

r/ = ^^2^ . (2.30) 
In order to evaluate the function h, one can set, using a conformal transformation, xq = oo, 

Xc — 0, Xqiq — 1 



X 



27a 

bo 



(1 - Xbo 



xixt,-xpr^^^x-^^^xf^\ (2.31) 



and therefore, 



Hv,v) = / d^x^d'xpix^ - x^)^>+2^-^(l - x^r^''^\v - Z;3r'^'"x-27.^-27. 

\ 27 / \ 27„ 



A I' 

bo J J 



d^x^d^xp ■■■ . (2.32) 

^XQa^XQc J \XaQcXQibo J J 

Remarkably, the same integral as in ( |2.15| ) appears, as a consequence of the symmetry 
between the limits xqi 00 and xqi found in [13]. First consider the case of large 
impact parameters b = |xoaol ~ \xobo\ ~ I^^OcI ^ r (see, fig. |^(b) and related calculations 
in [13,31]). Then rj is approximately a phase rj « e*^ where 6 is the relative angle as before. 
We find^ 

/ « 2c,f{i, la, la) (^) 9{e) , (2.34) 



^In fact, this result can be reac hed f rom ( 2.16 ) via a conformal transformation thanks to the conformal 
invariance of the original integral ( 2.13 ). Consider a SL(2,C) transformation 

(2.33) 



2-1/6 



Under this, one has xqi — > a^oi ~ xoi/a;oa;i « 4/2:01, Xc = ~* x'^ = b, Xag = r —* x'^^^ = 6/(1 — 6r), and 
Xaoc = r ^ r' ^ b^r. Therefore, 4r/xoi = x'q 
transformation does not change the angle 0. 



Xagc = r ^ r' ^ b^r. Therefore, 4r/xoi = Xq^t' /b^ as expected. Note finally that by definition a conformal 
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and 

r(2) ~ £ dy I djd^adjb {^y^ e^(7)y+(x(7.)+x(7.))(v--.) . (2.35) 

Again, the saddle points are given by 7a = 7fe = 1/2, and we have the pole at xil) = 
2x(l/2). On the other hand, the single scattering amplitude at large impact parameter is 

r(xoi,r,6)~a,(^)'%>^(^)^, (2.36) 

where 7 is the solution to 

X'(7)i" = ln^. (2.37) 

Taking the ratio, we find 

7^(2) / u2 \ 2(27-7) 
R=I11^ ±_] ,2(x(l/2)-x{7))l^ . (2.38) 

So in this case the correlation R decreases as either xqi or r (or both) is increased (keeping 
xoi,r < h). 

In order to exhibit a symmetry with respect to the large dipole case, let us look at the 
case of small impact parameters, typically, 6 ~ r ^ xqi- We find 

/~(^)^\ (2.39) 



while 



with 7 determined from 



r(xoi,r,6)~a,(^)'%>^(^)^, (2.40) 



y(^)y = ln^, (2.41) 

^01 



so that 



/ \ 2(27-7) 

— ) e2(x(l/2)-x(7))r_ (2.42) 

Compare with ( ^.27 ). As xqi increases, while keeping xqi <C r, the correlation decreases. 
From the limiting behaviors, (|2.27|) (xqi ^ r) and ( 2.42| ) (xqi ^ r), we see that R is 
enhanced when the asymmetry [xqi vs. r) is large, and it presumably takes a minimum 
value around xqi ~ r. 
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2.3 Estimates and comments 



Regarding the rapidity dependence, we note that 7 — > 1/2 as K — > 00. Thus for large 
Y, the coefficient multiplying Y in the exponent in ( ^.271) and ( |2.42| ) tends to zero. For 
a fixed Y, this coefficient again tends to zero when xqi ^ r, as can be seen from (|2^ ) 



and ( 2.41 ). Therefore the results ( p. 27 ) and ( 2.42 ) predict that the correlation R decreases 



faster with Y when XQi/r ^ 1 and XQi/r <C 1, while if we extrapolate our results towards 
the symmetric limit xqi ~ r, we see that R is almost constant in Y. 

From p36| ) we can guess that 7 is quite close to 1/2. Let us therefore set 7 = 1/2 + e 



and expand the BFKL eigenfunction to linear order in e. One then finds that 

^ :ln^= In^, (2.43) 



V'"(l/2)a,y r 14C(3)a,y 

where C(3) ^ 1.2. If xoi/r = 2 we then find, for = 0.2, e 0.21/Y, and thus for y = 8 
we have e ~ 0.03, while for y = 12 we find e ~ 0.02. For 2(27 — 7) we then find the 
values 0.46 and 0.43 for Y = 8 and 12 respectively. If instead XQi/r = 40 we find e « 0.14, 
2(27 - 7) « 0.91 and e « 0.09, 2(27 - ^ 0.73 for F = 8 and 12 respectively. For 
this values of 7 we also note that the exponent multiplying Y in ( p. 271) is quite small, for 
7 = 0.64 it is 0.14 while for 7 = 0.59 it is 0.06 (all these estimates are valid for = 0.2). 

Thus if, for a fixed Y, we try to fit as a function of XQi/r using a single effective 
power, uj, we would expect this fit to give a too strong increase close to the minimum, 
XQi/r ~ 1, whereas it should give a too slow increase further away from the minimum. As 
2(27 — 7) varies stronger for smaller Y, we would expect the fit to work better for higher 
Y. We would also expect lo to be larger for smaller Y. 

In the next section we will see that these analytical estimates are all in quite good 
agreement with the numerical results. In particular, the numerical analysis will confirm 
that the minimum of R (for zero impact parameter) occurs at xqi ~ r. Moreover, the 
estimates for 7 given above agree very well with the numerical results, and also the Y 
dependence turns out to be correct. 

Before moving on to the numerical analysis, we would like to address one more point. 
So far we have been able to make analytic estimates only for specific configurations. In 
particular, we assumed that the dipoles Xa^c and x^qc are more or less equal in size. In 
going from ( |1.2| ) to ( |1.1[ ), however, the question is whether the replacement 

dhMix,y,z)-TPix,z;z,y) ^ j cP z M{x,y, z) ■ Ty{x, z)TY{z,y) , (2.44) 

is valid. (We have here returned to the notation used in the introduction using x, y and z.) 
What we have shown above is that T'^^)(x, z] z, y) ^ T{x, z)T{z, y) for some specific regions 
of z, and also for specific relations between (x, y) and the target, but this is not sufficient 
to see the integrated effect of the correlation. Although one can use the MC code to do the 
integration over z, this can be quite time consuming. Leaving the numerical integration 
for future work, we here crudely identify the configurations which dominate the integral in 
( |1 . 1| ) . Consider the large parent case where |x — ?/| <S xqi and assume that |x — y| is smaller 
than the saturation length Q^^. This means that we may set T{x,y) = {x — y)'^Ql. (We 
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could also introduce an anomalous dimension 7 7^ 1 but this is not essential.) We then 
divide the integral into three regions: 

• Region A: \x — z|, \v — z\ <\x — y\ . 

• Region B: \x — y\ <\x — z\ \z — y\ < . 

• Region C: <\x — z\ \z — y\ . 
In region A we have 

d'^ ^,l\y^y'_,y [i^ - ^?Ql + - yfQl - - y?Ql - - ^?Ql ■ - v?Ql 

^{x-yfQl. (2.45) 

(Note that there is no logarithmic singularity at either z = x or z = y.) In region B we 
instead have 

' -{x - yfQl + 2{x - zfQl - (x - zfQ^ 



^{x-yf [ Sz-^—^2{x-zfQl 



while in region C we have 

|^d^.fc^|-(x - yfQl + l'^^{x-yfQl (2.47) 

where the integral is dominated by the lower limit |x — z| ~ l/Qs- Thus for a small 
projectile which has not yet reached saturation \x — y\ <^ Qj^^ the dominant contribution 
comes from region B where we indeed have \x — z\ \z — y\. As \x — y\ ^ 1/Qs, region 
B shrinks, and the dominant region is simply |a; — z| ~ |z — y| ^ \x — y\. Therefore, we 
expect that the configurations we are using are relevant, and the large correlation found 
there should survive after integrating over z in the evolution equation. 

3. Numerical Approach 

3.1 Outline of the approach 

In this section we will perform a numerical analysis to compute the quantities T^^^ and 
(r)^. This can be done rather easily in a Monte Carlo implementation of the dipole model, 
and we will here use the C-|— |- code developed in [18]. The calculation we will perform is 
straightforward, no matter which configuration we have. Recall that the definitions of T 
and T(2) are 

Tyi..y) = lM.M.Mu.v)nyiu,v). (3.1) 

T^\xi,yi;x2,y2)= / (fui(fvi (fu2(fv2Ao{xi,yi\ui,vi)Ao{x2,y2\u2,V2)n^Y\'^i^'"i'^'^2,V2) 
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where is the elementary dipole-dipole scattering amphtude. (The second term on the 



right hand side of (3.2) represents scattering of two dipoles off the same dipole in the 
target.) Starting from any initial dipole distribution, the MC code evolves the initial state 
up to a given value of Y, after which one can calculate all possible scatterings between the 



dipoles. The Monte Carlo estimate of equation ( |3.2| ) is simply given by 

TMci^i,yr,X2,y2) = j-j—Y^ ^ Ao{xi,yi\ui,Vi) ■ Ao{x2,y2\uj,Vj) , (3.3) 

n=i ijer„ 

where r„ is the configuration of the evolved target for the nth event. Writing X^t j ~ 
Z^i^j + X^i see that ( |3.3| ) contains both contributions in ( |3.2D . In writing this formula 
we only evolved the target but we can obviously do the computation in any given frame. 
Similarly the product T{xi,yi) ■ T{x2,y2) is calcuated as 



TMc{xi,yi) ■TMc{x2,y2) = Y Ao{xi,yi\ui,Vi) ■ — Ao{x2,y2\ui,Vi) 



n=l i&Fn 

(3.4) 

In the next section we will start by checking the predictions from [13] as stated in 



equations (1.4) and ( |2.7] ). As in the analytical approach we consider a target which initially 
consists of a single dipole {xq,xi) (for the numerical calculation we could start from any 
configuration if we so wish) For the configurations in [13], the phenomenologically more 
relevant configuration is the one in which the target xqi is much larger than the projectile 
dipoles. We fix the projectile dipoles to have the same size, r = Xa^ai = Xb^bi (for the above 
formulas this means we have xi = Xao,yi = Xa^,X2 = Xbo,y2 = a^fej, while the distance 
between them, Xab, will be varied. 

For the BK configurations, we have X2 = yi = Xc, and again we fix the two projectile 
dipoles to have the same size, r = \xao — Xc\ = \xc — Xbo\- The target dipole {xo,xi) is 
placed at zero impact parameter, as in figure ^ (a), while its orientation is chosen randomly 
for each event. We will always keep XaQ,XbQ and Xc fixed while we vary xqi and the impact 
parameter. 

One technical point is that one has to introduce a cutoff, p, for the minimal size of 
dipoles generated during the evolution since the dipole kernel A4{x, y, z) diverges ai z = x 
and z = y. Such a cutoff explicitly breaks conformal symmetry, and one should therefore 
ideally choose a cutoff which is much smaller than the relevant scales (the initial dipole 
sizes) involved in the process. On the other hand, simulations with too small values of p 
are very time-consuming. If one is studying symmetric collisions r xoi, then the choice 
p = O.Olr = O.Olxoi is good enough. Choosing an even smaller p in this case is not useful 
since one is then wasting a lot of time to generate many very small dipoles which do not 
interact and do not contribute much to the scattering amplitude. However, here we wish to 
study the correlation as we vary xqi, and then the choice of p is more subtle. For example, 
for a very asymmetric collision, say xqi ^ lOOr, p has to be much smaller than O.OIxqi so 
that we do not suppress important dipoles with size of order r. Besides, in the absence of 
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Figure 4: The numerical results for the configurations described in (1.4) at F = 10, and for target 
of size xqi — 20r (left plot) and xqi = 30r (right plot). The MC results are shown as squares while 
the power-like fits to the results are shown as dashed lines. 



saturation effects, smallness of p is also required for the frame-independence of T^'^\ hence 
that of R. As a compromise between these requirements (reducing simulation time and 
ensuring frame-independence) we shall choose p{xqi) = 0.05 r throughout. With this choice 
we confirmed that the results presented in what follows are reasonably frame-independent 
even up to the center-of-mass frame. 

3.2 Results 

As mentioned above we start by checking the results from [13]. The target will be fixed at 
the origin, with random orientation, and the projectile dipoles are placed symmetrically 
along the horizontal axis, one on the positive axis and the other on the negative axis, with 
random orientations. We choose as = 0.2 throughout, except in the running coupling case 
to be presented later. 

The results for this configuration are shown in figure |^. Here we choose xqi = 20 r in 
the left plot, and xqi = 30r in the right plot keeping xqi > Xab- The former case would 
in DIS correspond to a virtuality of ~ 60 GeV^. In both cases we also show fits of the 



form R = a/{xah + (^V ■ We thus confirm the power-like behavior in (1^), and also see 
that R converges to a finite value as Xab — > in agreement with the analytical prediction 
( |2.27D . For the left plot the fit gives the values j3 = 0.09 and 7 = 0.70 while for the right 



plot we get /? = 0.09 and 7 = 0.72. 

Next we turn to the BK configuration described in the previous section. In figure |5|, 
we plot i? as a function of XQi/r for y = 6, 8 and 10, at zero impact parameter. We can see 
a behavior of R consistent with the analytical formulas, equations ( |2.27| ) and ( |2.42| ). The 
minimum of R indeed occurs at xqi ~ r with the minimal value R ~ 1.5. For asymmetric 
configurations, R can easily reach values of order 10. Moreover, the powers extracted from 
figure ^ agree with the expectations from equations (2.21) and ( p. 26 ). For y = 6, a fit of 



the form {xQi/r)"^ gives the values a; = 0.55 in the region x^i/r = 1 — > 10, to = 0.93 in 
the region XQi/r = 6 — > 40, and uj = 1.22 in the region XQi/r = 40 — > 200. These values 
corresponds to 7 = 0.55, 0.64 and 0.72 respectively. If we instead calculate the power 
by calculating 7 using equation ( p. 43 ) at the points xoi/r = 5, 20 and 120 representing 
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Figure 5: The numerical results for R at Y = 6 (upper left plot), Y 
y = 10 (bottom plot) at zero impact parameter. 



8 (upper right plot) and 



the three regions above, we find the respective values 0.68, 0.96 and 1.32, in very good 
agreement with the numerical results. Similarly, for y = 8 we find the values uj = 0.52 and 
0.76 from fits in the first two regions above. This can be compared to the analytical result 
which gives 7 = 0.60 and 0.80. 

From our analysis in the previous section we know that R decreases as Y increases, 
and the rate of decrease is larger for asymmetric scattering. This tendency can be clearly 
observed, though the ratio R doggedly stays ^1.5. In the current simulation we cannot go 
to larger values of Y because the single dipole amplitude T for xoi ~ r reaches order unity 
around y = 10. Therefore, in the entire domain of Y values where our approach makes 
sense, the mean field approximation R = 1 is nowhere valid even in central collisions. Since 
this persists up to the onset of the strong scattering regime T ~ 0{1), it is unlikely that 
saturation effects immediately wash out the correlation. Rather, one has to carefully study 
the effect of correlations when solving nonlinear equations. 

Another, perhaps more striking consequence of the correlation emerging from our anal- 
ysis is that it makes the nonlinear term T^^) comparable to T even when T <C 1. For exam- 
ple, we have T = 0.023 for xqi = 40 r at y = 8, and in this case we see from figure ^ that 
R = 10. This means that T^^) = 0.0056, and thus T^^) /T = 0.24, so r(2) is not completely 
negligible as compared with T. For the more symmetric case xqi = 6r at y = 10 we have 
T = 0.39 while T^^) = 0.32 and R = 2.2, see again figure |. For xqi = 40r and y = 10 we 
instead have T = 0.059, while R = 6.9 and therefore T^"^^ /T = 0.41. Taken at face value. 
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Figure 6: The numerical results for R at nonzero impact parameter b 0, Y = 10 and xqi = lOr 
(left plot) and .tqi = 20r (right plot). 



these estimates suggests that one might have to include the nonlinear effects in the evolu- 
tion already in the dilute regime where T <C 1. We did not include such a back-reaction 
into our linear dipole evolution, and in this regard our analysis is not complete. This point 
certainly deserves further study. 

So far we have studied only configurations with zero impact parameter 6 = 0. At finite 
impact parameter the correlation becomes larger as suggested by ( 2.38| ). Of course if we 
think of xoi as representing the proton radius then one should be careful in interpreting 
results for b ^ xqi where confinement effects are certainly important. As a check of the 
analytical prediction, and also for the sake of demonstration, we nevertheless present 
some results when b > xqi. Figure ^ shows the b dependence of R for xoi/r = 10 and 
2^01/^ = 20. We see that R is almost constant as long as b is smaller than xqi and that it 
grows rapidly when b > xqi- 



Numerical simulation with a running coupling 



One of the non-leading effects which we can easily incorporate into the numerical 
simulation is the running coupling as has already been done in [18-20]. Although in this 
paper we mainly concentrate ourselves on the fixed coupling case, we would here like to 
briefly mention some of the preliminary results obtained when the running coupling is used. 

Technically, the inclusion of the running coupling is completely straightforward and 
we shall use the one- loop expression for Og, 



In^)ln(QVA^CD) 



(3.5) 



where we fix Aqcd = 0.22GeV. The running coupling enters both in the dipole evolution 
(as (Xs) and in the individual dipole-dipole scatterings (as a^). We will set A^'c = 3 and 
= 3 as in [19,20]. 

To avoid the IR singularity we shall freeze the coupling below a minimum scale Qmin 
corresponding to a maximum dipole size rmax = ^/Qmin- As in [20], we choose r^ax = 
3.5GeV~^. In [20], Oig was evaluated at the scale 1/Q = min(r, ri,r2) for the splitting 
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Figure 7: The crosses are the numerical estimates of R obtained using a running couphng at 
Y — 6, and for 6 = (left plot) and b = 5r (right plot). The squares are the corresponding fixed 
coupling results. 



r — > ri,r2, and this choice roughly follows from next -to-leading log (NLL) studies of the 
dipole evolution [32,33]. [See Section VII of [34] for a compact discussion.] Thus we continue 
to use this scale in the evolution of the dipole cascade. For the dipole-dipole interaction 
the correct choice of the scale is more subtle, and we here use the option described in [20]. 

In practice, simulations with the running coupling are quite time-consuming, and we 
have therefore not been able to check as many configurations as in the fixed coupling 
case. In figure we show the results obtained at y = 6 both at zero (left plot) and 
nonzero (right plot) impact parameter, together with the fixed coupling results. We see 
that R is somewhat reduced, but its minimum value is still around 1.5. We also see that 
the qualitative behavior of R does not change, the minimum again occurs when xqi ~ r 
although it is of course difficult to determine the exact behavior of R since we do not have 
enough data points. At y = 8 for 6 = 0, we find the value R = 1.5 at xqi = 2r, while in 
the fixed coupling case we found R = 1.6. For 6 = 5r, R reduces from 11.6 in the fixed 
coupling case to 9.4 in the running coupling case for the same configuration. 



4. Conclusions 

In this paper we have studied both analytically and numerically the correlations induced 
by the leading order BFKL dynamics in the high energy evolution of a dilute system (such 
as a proton). Our main analytical results are given in equations ( ^^ , (|2!38|) and (l2l^ ). 



All these results indicate that one should expect power-like correlations which lead to a 
strong violation of the factorization T^^^ ~ T • T. The analytical estimates have been 
demonstrated to be qualitatively correct by a numerical analysis with which we have also 
been able to quantitatively study the behavior of the ratio R = T^^^/T^. We have found 
that R is always larger than ~ 1.5 and it can easily reach ~ 0(10) when the asymmetry is 
large. 

Physical consequences of the correlation remain to be explored. The first and obvious 
intuition is that it opens an intriguing possibility of the 'grey disc' limit in which a scattering 
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amplitude saturates to a value less than 1.^ 



r^-<i. (4.1) 

However, since R is not a constant, and the nonlinear equations involve an integration over 
the transverse plane with a nontrivial weight, a more detailed analysis would be required 
in order to draw any conclusions. 

Another interesting problem is the interplay with the gluon number fluctuation which 
has attracted considerable attention lately (see [36] and references therein), but which 
has so far mostly been studied in simple toy models where the transverse dimensions are 
suppressed. Though it typically requires unrealistically large energies to see the impact of 
the gluon number fluctuation on the nonlinear evolution of large nuclei, this is probably 
not the case for a dilute target. The BFKL evolution generates a very strong number 
fluctuation as well as the transverse correlation in the dilute regime, and they can both 
affect the subsequent nonlinear evolution in significant ways. 

There is plenty of room for improvements in the Monte Carlo simulation itself. In order 
to make a quantitative prediction for realistic experiments, one should include various NLL 
corrections and saturation effects into the target evolution. They have been incorporated in 
the dipole model in [18-20]. Among them, we have in this paper included some results with 
the running coupling effect. Since our simulations have been limited in size, it is difficult 
to determine the exact behavior of R. What we have clearly observed, however, is that R 
is somewhat reduced from the fixed coupling case, but is still large. This suggests that the 
large correlation may not be totally attributed to conformal symmetry of the leading order 
BFKL, but rather is a robust feature of the QCD evolution in the linear regime. 

As mentioned in the introduction we would expect even larger correlations in the 
multiple scattering amplitudes T^^^ (p > 3) which enter the Balitsky hierarchy. In the 
dipole model, these amplitudes are directly related to the corresponding multiple dipole 
distributions n^^^ [23,24], but analytical results for them are scarce [37]. The numerical 
evaluation of these amplitudes is straightforward, although the calculation of T^^^ for large 
p would be time-consuming due to the need of good statistics. 
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